C Kolekce procedur na reseni soustav linearnich rovnic s pozitivne
C definitni matici Choleskeho rozkladem. Od J.Humlicka.
C
C  Tato verze pracuje s DVOJNASOBNOU PRESNOSTI.
C
C 1996-12-23 F.Hroch


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC      
      
      SUBROUTINE DLL(NMAX,N,A,D)
C PROCEDURA NA ROZKLAD REALNE SYMETRICKE POZITIVNE DEFINITNI MATICE
C NA DVE TROJUHELNIKOVE MATICE, PROMENNE JSOU V DOUBLE PRESNOSTI
C   'NMAX' JE DEKLAROVANA DIMENZE ROZKLADANE MATICE, 'N' SKUTECNA DIMENZE
C   'A' JE ROZKLADANA MATICE KTERA JE V HLAVNIM PROGRAMU DEKLAROVANA
C   JAKO A(NMAX+1,NMAX), TEDY MA O JEDEN RADEK VIC NEZ SLOUPEC,
C   ROZKLAD PAK PONECHA HORNI TROJUHELNIK A ROZLOZENA MATICE JE 
C   VE SPODNIM ALE POSUNUTA O RADEK DOLU TEDY MA PRVKY A_{I,J} = A(I+1,J),
C   NA MISTO DIAGONALNICH PRVKU JSOU ULOZENY JEJICH PREVRACENE HODNOTY 
C   A_{I,I} = 1.0/A(I+1,I)
C   PRI VOLANI STACI NADEFINOVAT POUZE HORNI TROJUHELNIK
C   'D' JE INDIKACE POZITIVNI DEFINITNOSTI, PRO POZITIVNE DEFINITNI
C   MATICI JE D=1.0 JINAK D=0.0, PROCEDURU JE MOZNE MODIFIKOVAT TAK,
C   ABY V PROMENNE 'D' BYL DETERMINANT, ALE PRI JEHO VYPOCTU MUZE DOJIT 
C   K PRETECENI
C      
      INTEGER N,NMAX,I,J,K
      DOUBLE PRECISION D,A(NMAX+1,NMAX)
      DOUBLE PRECISION S,DD
                       
      D =0.0
      DO I=1,N
        DO J=I,N
          S = A(I,J)
          DO K=1,I-1
            S = S - A(I+1,K)*A(J+1,K)
          ENDDO
          IF(I.EQ.J)THEN
            IF(S.LE.0D0) RETURN
            DD = 1D0/DSQRT(S)
            A(J+1,I) = DD
          ELSE
            A(J+1,I) = S*DD
          ENDIF
        ENDDO
      ENDDO
      D = 1.0
      END
                  
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE DLLSOL(NMAX,N,A,B,X)
C PROCEDURA NA RESENI SOUSTAVY LINEARNICH ROVNIC S MATICI ROZLOZENOU      
C PROCEDUROU LL, PROMENNE JSOU V DOUBLE PRESNOSTI
C  PROMENNE 'NMAX','N','A' MAJI STEJNY VYZNAM JAKO V LL
C  'B' JE VEKTOR PRAVYCH STRAN A NEMENI SE
C  'X' JE VYPOCTENE RESENI

      INTEGER N,NMAX,I,K
      DOUBLE PRECISION A(NMAX+1,NMAX),B(NMAX),X(NMAX)
      DOUBLE PRECISION S
      
      DO I=1,N
        S = B(I)
        DO K=1,I-1
          S = S-A(I+1,K)*X(K)
        ENDDO
        X(I) = S*A(I+1,I)  
      ENDDO
      DO I=N,1,-1
        S = X(I)
        DO K=I+1,N
          S = S-A(K+1,I)*X(K)
        ENDDO
        X(I) = S*A(I+1,I)    
      ENDDO
      END
        
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      SUBROUTINE DLLINV(NMAX,N,A)
C VYPOCET INVERZNI MATICE PRO REALNOU POZITIVNE DEFINITNI MATICI     
C PREDTIM ROZLOZENOU POMOCI PROCEDURY LL, PROMENNE JSOU V DOUBLE PRESNOSTI
C PROMENNE MAJI STEJNY VYZNAM JAKO V LL
C VYSLEDNA INVERTOVANA MATICE JE ULOZENA JEN JAKO DOLNI TROJUHELNIK
C A^{-1}_{I,J} = A(I+1,J), HORNI TROJUHELNIK SE NEMENI
      INTEGER N,NMAX,I,J,K
      DOUBLE PRECISION A(NMAX+1,NMAX)
      DOUBLE PRECISION S
      
      DO I=1,N
        DO J=I+1,N
          S = 0D0
          DO K=I,J-1
            S = S-A(J+1,K)*A(K+1,I)
          ENDDO
          A(J+1,I) = S*A(J+1,J)
        ENDDO
      ENDDO       
      DO I=1,N
        DO J=I,N
          S = 0D0
          DO K=J+1,N+1
            S = S+A(K,J)*A(K,I)
          ENDDO
          A(J+1,I) = S
        ENDDO
      ENDDO
      END
                  
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

